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I .  INTRODUCTION 


A  burst  of  debris  fragments  is  produced  when  an  AP  round  or 
shaped-charge  jet  perforates  the  armor  of  a  tank.  These  fragments 
can  damage  interior  components  and  injure  crew  members  and  thus 
contribute  to  the  disruption  of  the  proper  operation  of  the  vehicle. 
Consequently,  in  Army  vehicle  vulnerability  studies,  this  spall- 
burst  damage  must  be  evaluated  and  added  to  that  caused  by  the 
primary  round  or  jet  to  obtain  the  total  incapacitation  of  the 
vehicle . 

Experimental  spall-burst  data,  gathered  at  the  Ballistic 

12  3 

Research  Laboratory  (BRL)  ’  *  of  the  US  Army  Armament  Research  and 
Development  Command  (USARRADCOM) ,  were  given  to  Systems,  Science,  and 
3 

Software  (S  )  to  analyze  and  develop  into  a  spall-burst  model. 

4 

Their  primary  conclusion  was  that  the  available  data  were  so 
sparse  and  possessed  such  huge  scatter  that  the  final  development 
of  a  predictive  model  was  impossible  at  that  time.  However,  they 
did  conclude  that  the  distribution  of  fragments,  averaged  over 
several  similar  firings,  appeared  to  be  symmetrical  about  an  axis 
located  midway  between  the  normal  to  the  inner-armor  surface  and 
the  direction  of  the  perforating  round  or  jet  (Figure  1). 

Calculational  studies  have  been  conducted  at  BRL  to  determine 
the  increase  in  the  kill  rate  of  vehicles  when  the  contribution  of 


1 

B.  Rodin ,  "Documentatvon  of  Spall  Produced  by  a  3.2  Inch  Precision 
Shaped  Charge , "  BRL  Interim  Memorandum  Report  No.  226,  May  1974. 

2 

R.J.  Pipino,  R.C.  Brown ,  and  L.T.  Brown,  "Compilation  of  Data  on 
Spall  Produced  from  Armor  Plate  by  Detonation  of  a  Shaped  Charge," 
THOR  Technical  Note  No.  138,  June  1968. 

$ R.J .  Pipino  and  L.T.  Brown,  "Spall  Decks  and  Tapes-Identification 
and  Contents,  "  THOR/Falcon  Letter  Report  dated  November  9,  1971. 

4 

R.T.  Sedgwick,  P.L.  Anderson,  M.S.  Chawla,  C.R.  Hastings,  and 
L.J.  Walsh,  "Characterization  of  Behind-the-Armor  Debris  Resulting 
from  Shaped-Charge  Jet  Formation , "  BRL  Contract  Report  No.  249, 
July  1975. 


the  spall-burst  damage  is  included.  A.  !laferJ  has  calculated  the 

kill  rate  of  a  tank  exposed  to  shaped-charge  land  mines.  T.  Hafer6 
has  calculated  the  kill  rate  of  a  tank  exposed  to  shaped-charge 
rounds.  Each  of  these  calculations  was  performed  using  the  VAST7 
assemblage  of  computer  programs.  A  significant  contribution  by 
spall  bursts  to  the  incapacitation  of  vehicles  was  reported  by  each 
study. 

The  components  which  are  essential  to  a  specified  operation  of 
the  vehicle  and  are  susceptible  to  being  disabled  by  the  direct 
impact  of  representative  projectiles  are  identified  as  critical. 

The  remaining  components  of  the  vehicle  are  identified  as  inert 
and  participate  in  vulnerability  studies  only  insofar  as  they  may 
shield  components.  The  fractional  loss  of  capability  (incapacitation) 
by  a  critical  component  is  given  by  the  incapacitation  version  of  the 

g 

disablement  function.  An  incapacitation  function  can  often  be 
approximated  by  a  kill  function  which  gives  the  probability  of  a 

8  9 

total  loss  of  capability  (kill)  per  impact.  ’  This  study  uses 
only  the  kill  version  of  the  disablement  function. 


l/1.  Hafer ,  "Analytical  Evaluation  of  Spall  Suppression  in  the  Soviet 
T-55  Tank  as  a  Countermeasure  to  the  U.S.  XM-70,"  Mine-Countermeasure 
Quarterly  Symposium ,  Picatinny  Arsenal,  New  Jersey,  1  October  1975. 

'  T.  Hafer,  "The  SHUTE  Internal  Point  Burst  Vulnerability  Assessment 
Program,"  Proceedings  of  the  2nd  Symposium  on  Vulnerability  and 
Survivability,  Volume  II,  January  1977. 

"Vulnerable  Areas  for  Surface  Targets  (VAST) ;  An  Internal,  Point- 
Burst  Model,"  compiled  by  C.  Nail,  Computer  Sciences  Corporation, 
Contract  DAAK40-78-D-003/P0009,  September  1978. 

Q 

W.B.  Beverly,  "The  Disablement  Functions  of  a  Critical  Component 
and  their  Use  in  Vehicle  Ballistic  Vulnerability  Calculations, "  to 
be  published  as  a  USAARADCOM/BRL  Report. 

g 

L.  Kruse  and  P.  Brizzolara ,  "An  Analytical  Method  for  Deriving 
Conditional  Probabilities  of  Kill  for  Target  Components ,"  US  Army 
Ballistic  Research  Laboratories  Report  No.  1563,  December  1971. 
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In  the  stochastic  model  of  Beverly,  the  descriptions  (states)  of 
representative  samples  of  fragments  which  have  perforated  a  given  amount 
of  material  are  assumed  to  he  given  by  probability  density  functions 
(I’D!-)  of  incomplete  sets  of  parameters.  The  killing  probability  of  a 
spall  burst  (Figure  2)  described  by  B*(Sq)  is  calculated  in  terms  of 

the  expectation  of  the  POP  of  each  impacting  fragment  and  the  kill 
function  R*(;tj)  of  the  component.  In  analytic  form,  the  kill  prob¬ 
ability  I’  of  the  spall  burst  is  given  by 

P  =  1  -  e'x  (1A) 


where  \  is  given  by  either 

x  ^LbVoItVo.q  -PViKo)] 

P  (oi|a0)R  (o1)dA1  dA0 ,  (IB) 


or  the  more  concise  form 


The  quantities  used  in  the  preceding  equations  are  defined  as: 

a  =  a  set  of  parameters  which  incompletely  describes 
a  source  fragment, 

aj  =  a  set  of  parameters  which  incompletely  describes 
a  residual  fragment  at  impact  with  a  critical 
component , 

-*•  =  a  set  of  parameters  which  describes  the  medium 

®  between  the  fragment  origin  and  the  critical 

component , 

1°W.B.  Beverly ,  "The  Stochastic  Representation  of  the  Perforation  of 
Barriers  by  Steel  Fragments  in  Military  Vehicle  Ballistic  Vulner¬ 
ability  Studies to  be  published  as  a  USAA RRADCOM/BRL  Report. 


10 


R* (a j )  =  the  average  killing  probability  o£  a  representative 

sample  of  fragments  described  by  a^,  as  they  impact 
the  critical  component, 

-V  — ► 

P*(a i/aQ)  =  The  PDF  of  residual  fragment  states  derived  from  a 
source  fragment  described  by  aQ, 

T* [a  ,g-^P*(a^/a  )]  =  The  operation  which  describes  the  average 

transport  of  a  representative  sample  of 
source  fragments ,  all  described  by  to  a 

possible  impact  with  the  critical  component. 

The  spall-burst  source  term  normalizes  to  the  number  of  signifi¬ 
cant  fragments  N;  i.e.. 


_^B*(a^)dAos  N.  (ID) 

The  incomplete  parameter  sets  aQ  and  a^  are  composed  of 

°0  =  [r0  *  Po  Q0  ^7  *  *  *  *  (  ao  V  ]  »  (IE) 

al  =  [ri»Pl?(al)7»*  ]>  (IF) 

where  rQ  and  r^  are  position  vectors  which  locate  the  fragments,  Pq 

and  Pj  are  the  directions  of  the  fragment,  (aQ)^  and  (a^)y  are  the 

seventh  components  of  a^  and  a^,  and  K  is  the  dimensionality  of  both 

a  and  a  j .  We  will  use  the  more  explicit  form  of  equation  IB  in  the 

following  discussion  since  this  form  can  be  more  easily  related  to 
the  Monte  Carlo  solution  which  is  given  later  in  this  report. 
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The  integral  in  equation  IB  is  identified  here  as  the  kill 
integral.  The  star  superscript,  sometimes  used  to  identify  adjoint 

quantities,  is  used  here  to  identify  problem-dependent  quantities. 

The  dagger  superscript  is  introduced  later  in  this  study  to  identify 
the  adjoint  transport  operation. 

We  have  assumed  in  equation  IB  that  the  shape  and  orientation 
of  each  fragment  are  described  in  some  detail.  Since  this  detail  can 
not  be  predicted  by  the  current  spall-burst  and  fragment-transport 
models,  we  must  modify  equation  IB  to  form  which  is  compatible  with 
the  present  level  of  detail  in  these  models.  The  nomenclature  of 
equation  IB  is  retained  insofar  as  possible,  but  the  form  of  the 
functions  and  operations  will  change  as  their  parameter  lists  are 
shortened.  The  transport  operation  is  replaced  by  a  simpler  operation 
which  merely  predicts  the  average  mass  and  velocity  of  the  largest 
residual  fragment.  The  spall-burst  source  term  and  the  critical 
component  kill  function  are  averaged  over  a  representative  sample  of 
shapes  and  orientations  and  are  also  given  in  units  of  the  mass  and 
velocity  of  the  fragment. 

We  will  apply  four  additional  approximations  to  the  description 
of  spall-burst  phenomena  to  bring  vulnerability  calculations  to  a 
reasonable  level  of  effort.  The  effect  of  these  approximations  on 
the  accuracy  of  vulnerability  predictions  is  not  investigated  here 
but  should  be  evaluated  in  future  studies.  The  first  and  third  of 
the  approximations  are  suggested  in  Reference  4.  These  approximations 
are  briefly  tabulated  as: 

1.  All  fragments  in  a  burst  (source  fragments)  originate  at 
the  same  point. 

2.  All  fragments  follow  straight  paths  from  their  birth  to 
their  interaction  with  the  critical  component. 

3.  The  probability  density  functions  (PDF)  describing  an 
average  burst  are  symmetrical  about  the  axis  described 
in  Figure  1. 

4.  The  PDF  of  states  for  a  fragment  emerging  from  a  transport 
operation  is  replaced  by  the  mean  value  of  the  used  parameters. 
The  uncertainty  introduced  into  vulnerability  calculations  by 

this  approximation  is  being  investigated  at  BRL.^ 


22 

W.B.  Beverly ,  "The  Use  of  Correlated  Sampling  in  the  Monte  Carlo 
Calculation  of  Ballistic -Vulnerability  Predictions  for  Military 
Vehicles , "  to  be  published  as  a  USAARRADCOM/BRL  Technical  Report. 
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RT  COMPONENT 


Applying  these  approximations,  equation  IB  becomes 

*00  *00 

A  •/  /  /  B*<ro>mo»  vo>®> 

*0  ^0  •'o 

T*(^,mo,v8,0,g  —7,,^,,^,) 
R*{f’,m|1vl)dm0dv0de ,  (g) 


where 

mQ  =  the  mass  of  a  source  fragment, 

vq  =  the  velocity  of  a  source  fragment, 

9  =  the  angle  made  by  the  fragment  path  with  the  spall-burst 

axis  of  symmetry, 

nij  =  the  average  expected  mass  of  the  residual  fragment  at 

impact  with  a  critical  component  for  all  source  fragments 
having  mass  mp  and  velocity  v  , 

Vj  =  the  average  expected  velocity  of  the  residual  fragment 
at  impact  with  a  critical  component  for  all  source 
fragments  having  mass  mQ  and  velocity  v  . 

The  transport  of  steel  fragments  through  inert  material 
(non-critical  components)  is  currently  described  by  the  THOR 
12 

parametric  equations  given  below.  The  coefficients  used  in 
these  equations  have  been  evaluated  for  different  materials  using 
a  modified  least-squares  technique  to  fit  the  data  gathered  from 
firing  experiments.  The  coefficients  are  averaged  over  fragment 


12 

"The  Resistance  of  Various  Metallic  Materials  to  Perforation  by 
Steel  Fragments ;  Empirical  Relationships  for  Fragment  Residual 
Velocity  and  Resident  Weight,"  Project  THOR  Technical  Report 
No.  47,  April  1961,  The  Johns  Hopkins  University,  Baltimore, 
Maryland. 
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shapes  to  the  extent  that  separate  sets  of  coefficients  are  calcu¬ 
lated  for  chunky  fragments  and  for  fragments  having  no  particular 
shape.  The  THOR  perforation  equations  are: 

*  C6  C7  C8  C9  CO  A  % 

m,=  m0-IO  (AT)  m0(Sec*7)  v0,  (3A) 

Cl  C2  C3  C4  C5 

v,=ve-IO  (AT)  m0(Seci,)  v0,  (3B) 

where 

=  the  average  mass  (grains)  of  the  largest  perforating 
fragment, 

v  =  the  average  velocity  (ft/sec)  of  the  residual  fragment, 

m  =  the  mass  of  a  source  projectile  (the  projectiles  used 
0  here  are  not  actual  fragments.), 

vq  =  the  velocity  of  the  source  projectile. 

Cl,'"", CO  =  the  fitting  coefficients  (A  set  of  values  for  the 

coefficients  was  derived  for  each  target  material.), 

A  =  average  impact  area  of  the  fragment  (in2)  , 

T  =  target  thickness  (in), 

n  =  the  angle  between  the  fragment  shot  line  and  the 
normal  to  the  target  material. 

The  general  outlines  of  three  different  Monte  Carlo  methods 
for  evaluating  the  kill  integral  in  equation  2  are  developed  in 
Section  II  of  this  report.  The  necessary  transformations  and 
changes  are  developed  and  explained  for  each  method.  The  three 
methods  are: 

1.  The  Forward  Transport  of  Sample  Fragments  from  the  Spall- 
Burst  Origin  to  A  Possible  Impact  with  the  Critical  Component. 
Importance  sampling  (Appendix  A)  is  not  used  in  the  picking  of  sample 
path  rays.  The  PDF  of  sample  fragments  in  this  method  is  approxi¬ 
mately  that  of  fragments  in  an  actual  burst. 

2.  The  Forward  Transport  of  Sample  Fragments  from  the  Spall- 
Burst  Origin  to  A  Possible  Impact  with  the  Critical  Component. 

13 

A  method  of  importance  sampling  was  suggested  by  Reisinger  where 


M.  Reisinger ,  USAERADCOM/BRL ,  Private  Communication. 
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the  path  of  a  sample  fragment  is  constructed  through  a  point  picked 
with  equal  probability  at  any  location  within  a  rectangular  volume 
which  barely  encloses  the  critical  component.  An  approximation  im¬ 
plementation  of  this  method  is  used  in  the  current  BRL  vehicle 

14 

vulnerability  methodology.  An  exact  analytical  representation  of 
this  method  is  derived  here  and  its  use  in  vulnerability  studies  is 
outlined.  The  PDF  of  sample  fragments  in  this  method  differs  from 
the  PDF  of  fragments  in  an  actual  burst. 

3.  The  Adjoint  Transport  of  Sample  Residual  Fragments  from 
the  Outer  Surface  of  the  Critical  Component  to  the  Spall -Burst 
Origin.  The  importance  sampling  used  in  the  second  method  is  also 
used  here.  The  sample  events  used  here  are  calculational  in  their 
nature  and  do  not  correspond  to  any  actual  events. 

In  Section  III  of  this  report,  a  test  problem  is  devised  and 
its  Monte  Carlo  solution  is  calculated  using  each  of  the  preceding 
approaches.  In  each  solution,  the  transport  of  fragments  is  assumed 
to  be  described  by  the  THOR  equations  3A  and  3B.  In  these  solutions, 
the  geometrical  and  distributional  symmetries  of  the  test  problem  are 
utilized  to  simplify  the  calculations.  However,  the  solutions  demon¬ 
strate  the  essential  features  of  vehicle  vulnerability  calculations 
for  spall  bursts. 


II.  THE  MONTE  CARLO  CALCULATION  OF  THE  KILL  PROBABILITY 
OF  A  CRITICAL  COMPONENT  BY  A  BURST  OF  FRAGMENTS 

A.  The  Forward  Transport  of  Fragments  Without  Importance  Sampling 


The  general  solution  of  equations  2A  and  2B  are  outlined  below. 
In  this  solution,  each  calculational  event  will  correlate,  within 
the  accuracy  of  the  vehicle  vulnerability  model, *0  to  a  real  event. 

1.  Pick  a  set  of  parameter  values  for  a  source  fragment  by 

sampling  the  source  term  (Appendix  A).  These  parameter  values  are 

identified  as  fr  , (m  ) . (v  ).,0.1  where  i  is  the  running  index  over 
1  o  o  i»  o  1  iJ  ^  & 

sample  fragments.  The  burst  origin  r^  is  included  as  an  argument 

although  all  fragments  in  a  burst  are  assumed  to  originate  at  the 
same  point. 

2.  Transport  the  fragment  through  any  intervening  inert 
components  to  a  possible  impact  with  the  critical  component.  The 


L.  Losie  and  T.  Ha  fen,  "HIP”,  to  be  published  as  a  USAARRADCOM/ 
BRL  Report. 
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parameter  values  of  an  impacting  fragment  are  identified  as 

3.  Calculate  the  score  for  the  sample  event  as 

p*  t(rp  (m^)  . ,  (v-^)  .  ] ,  In  this  instance,  the  score  is  the  kill 

probability  of  the  impact.  A  fragment  which  does  not  impact  the 
component  is  given  a  score  of  0. 

4.  Repeat  steps  1  through  3  for  a  total  of  I  sample  fragments. 
The  value  of  the  kill  integral  is  then  estimated  using 


X  *  (N  2  Si)  /  I  *  N  S. 


(4A) 


A  bar  is  placed  over  a  variable  to  identify  a  mean  value.  We  have 
assumed  in  equation  4A  that  a  Monte  Carlo  estimate  is  made  of  N  in 
a  secondary  calculation  which  parallels  the  main  calculation.  The 
details  of  such  a  secondary  calculation  are  not  discussed  in  this 
report  but  can  be  easily  derived  by  using  the  procedures  outlined 
in  Appendix  A. 

5.  The  quantity  X  will  converge  toward  the  true  value  of  X  as  an 
increasing  number  of  fragment  histories  are  conducted  and  their  kill 
probabilities  are  calculated  and  assimilated  into  equation  4A.  The 
standard  deviation  of  S,  used  as  a  measure  of  the  statistical  accuracy 

of  S,  is  given  by15 


It  —  t  I/- 

»5-[I*-15  ]'  , 

1  Hi-i)  1 

and  the  standard  deviation  of  X  is  given  by 

8  X  -[(3  8  N)*+(N  8  S)  ]  , 


(4B) 


(40 


o 

Y.  Beers,  "Introduction  to  the  Theory  of  Errors,"  Addison-Wesley 
Publishing  Company ,  Inc.,  1952. 


I, 


where 


6N  *  the  standard  deviation  of  N", 

6S  =  the  standard  deviation  of  !?, 

6 X  *  the  standard  deviation  of  X. 

b .  An  estimate  of  the  kill  probability  of  the  component  by  the 
burst  of  fragments  is  given  by 

PK=l-e'x.  (IA) 

The  standard  deviation  of  the  kill-probability  estimate  is  given  by 

8PK=e'X8\  ,  (5) 

where  6P  is  the  standard  deviation  of  F. 


Step  6  completes  a  brief  outline  of  the  calculation  of  an  estimate 
of  P  where  sample  fragments  are  picked  by  sampling  the  unbiased  spall- 
burst  term.  In  this  solution,  each  calculated  score  S  is  an  approxi¬ 
mation  of  the  average  kill  probability  of  all  fragment*  described  by 


B.  The  Forward  Transport  of  Fragments  With  Importance  Sampling 


The  preceding  procedure  will  often  result  in  the  picking  of  a 
large  number  of  sample-fragment  paths  which  miss  the  critical  component. 
In  these  cases,  the  calculations!  effort  wasted  by  the  generation  and 
tracking  of  misses  can  often  be  reduced  by  changing  the  direction 
variables  of  source  fragments  to  new  variables  whose  region  of  inte¬ 
gration  permits  the  picking  of  sample  fragments  which  impact  the 
critical  component  with  a  higher  probability  (importance  sampling). 

We  choose  these  new  variables,  whose  region  of  integration  RP  is  the 
interior  of  the  rectangular  parallelepiped  which  barely  encloses  the 
critica^  component  (Figure  3),  to  be  the  coordinates  of  the  position 
vector  r7. 

A  sample -fragment  path  is  constructed  by  first  picking  a  point  r2 
with  equal  probability  at  any  location  in  the  region  RP.  A  ray  is 
then  constructed  from  the  burst  origin  t  through  the  sample  point  ?2- 
An  alternate  region  of  integration  RC  can  be  defined  as  that  part  of 
region  RP  which  will  always  yield  rays  which  intersect  the  critical 
component  (Figure  4).  The  quantities  used  in  the  transformation  are 
illustrated  in  Figure  5  and  are  developed  below  in  Equations  6A-6H: 


IS 


RPP  WALL  &  FRAGMENT  PATH  INTERSECTION 


Figure  3.  An  Example  of  a  Critical  Component  Inside  An  RPP 


The  new  quantities  introduced  in  the  preceding  equations  are  defined 

as : 

r  =  the  distance  along  a  ray  from  the  spall-burst  origin  to  the 
first  intersection  with  the  RPP, 

R  =  the  length  of  the  ray  intercepted  by  the  RPP, 

dS  *  an  infinitesimal  area  on  the  sphere  whose  center  is  the  spall- 
burst  origin  and  whose  radius  is  r  (This  S  is  used  only  in 
equations  6C  and  6G  to  represent  area  and  should  not  be 
confused  with  the  S  used  elsewhere  to  represent  scores.)* 

d£I  =  the  solid  angle  subtended  by  dS  at  the  spall-burst  origin, 

r 2  =  the  position  vector  which  locates  a  point  in  region  RP  or  RC* 

dV  =  the  infinitesimal  volume  in  the  RPP  which  is  associated  with 

dS  (Figure  5), 

Sc  =  the  surface  of  the  critical  component  which  can  be  impacted 
by  fragments  from  the  spall  burst, 

RP  =  the  region  enclosed  by  the  RPP, 

RC  =  the  region  inside  the  RPP  which  barely  encloses  all  possible 
fragment  rays  which  intersect  the  critical  component, 

VRp  =  the  volume  of  the  region  RP, 

VR(,  =  the  volume  of  the  region  RC. 

*  4 

In  its  present  form,  the  source  term  N  [ro,mo,vo, 0(r2)]  *n 

equation  6E  is  not  tractable  to  the  efficient  picking  of  sample  frag¬ 
ments.  Using  the  procedure  outlined  in  Appendix  A,  we  will  first 
change  the  kill  integral  to  a  form  where  N*  is  given  in  terms  of  the 
one  dimensional  function  G(r,)  and  the  two  dimensional  function 

Q(VWr2  Sample  paths  can  now  be  picked  by  sampling  1/VRC , 

and  their  associated  mass  and  velocity  values  are  picked  by  sampling  Q. 
These  changes  in  the  kill  integral  are  described  by: 


23 


-g%{r  /7aQ(mP.<o.^^)] 

rcjrc  r2F(r, R)  l«rc  *'cro 

T*[rp,nn0  ,v0  ,fl(r2),  9  —  7j(  ^ ) ,  m, ,  v] 

R*[?(^)i™ii5]dm0dv0dV2’|  dV2 

^VrC  ' 

_  ,®,®  -  — 

G[7~2 )  J  N*[l^ ,  mQ  ( vQ  ,0( r2  )]dmQ dv0  , 

Q  (  mo  1  v0  ’  r2  1  r2  ) 

-CD  -CD  ^  __  ' 

J0  J0N*[~0'm 0,v0'^(?2)]clm0c,v0 

S(r2“r2'):  ^[(r2  )|-(r2' )j]  S[(r2)2-(r2)2]B[(r2)3-(r2)3]  i 

where 

r7  =  a  dummy  value  for  r2, 
dV^  =  a  dummy  value  for  dV2> 

RC'  =  a  dummy  representation  of  region  RC, 


(7A) 

(7B) 


(7C) 


(7D) 
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6 [ (r-)  ,  - (ri)  ]  =  the  Dirac  delta  function  for  the  x-component 
1  1  of  r2  and  r£. 

A  brief  outline  of  a  Monte  Carlo  evaluation  of  the  preceding  kill 
integral  is  given  below: 

1.  Construct  an  RPP  which  barely  encloses  the  critical  component. 

2.  Pick  a  point  with  equal  probability  at  any  location  within  the 
enclosing  RPP.  A  running  tally  U  is  maintained  of  the  sample  points, 
and  is  used  at  the  end  of  all  sampling  to  calculate  an  estimate  of  VR(.. 

3.  Construct  a  ray  from  the  spall -burst  origin  through  the  sample 
point  and  determine  if  the  ray  intersects  the  critical  component. 
Discard  any  ray  which  misses  the  component  and  pick  another  ray  until 
an  intersection  is  obtained.  A  sample  point  which  leads  to  an  inter¬ 
section  is  identified  as  (?' ) ^  where  i  is  the  sample- fragment  index. 

4.  Calculate  the  quantities  r. ,  R. ,  (Figure  5)  and  0.  for  the 
sample  ray. 

5.  Pick  values  for  the  mass  and  velocity  of  the  sample  fragment 

by  sampling  Q[mo>vo,r2,  (r')  i] .  In  practice,  Q  is  changed  to  one¬ 

dimensional  functions  in  a  manner  similar  to  that  used  in  the  picking 
of  path  rays,  but  we  will  not  outline  the  operation  here.  The  full 
set  of  values  is  identified  as  [ , CmQ ) ^  » C vq ) i , © i ] . 

6.  Calculate  the  transport  of  the  sample  fragment  from  the  spall- 
burst  origin  to  a  possible  impact  with  the  critical  component.  The 
parameters  of  the  residual  fragment  at  impact  are  identified  as 

tcVi'C^i’Cv^il- 

7.  Calculate  the  score  of  the  event.  A  fragment  which  is 
defeated  before  impacting  the  component  is  given  a  score  of  0.  The 
score  of  a  fragment  which  impacts  the  component  is  given  by 


R*{n[(^)iMmi]i.h]i}G[(^)i] 

r?F[ri,R|] 


(8A) 


8.  Repeat  steps  2-7  until  a  total  of  I  sample  fragments  have  been 
picked  whose  paths  intersect  the  critical  component. 


*  * 


9.  Calculate  an  estimate  of  X  using 

(8B) 

V  =  ( i  V)  /u, 

(SC) 

where  U  is  the  number  of  sample  points  picked  in  region  RP. 

10.  Using  the  appropriate  form  of  equation  4B,  calculate_the 
standard  deviation  of  VR(_.  and  S.  The  standard  deviation  of  X  is 
then  calculated  using 

»x.[^»sf  +  ,sss«fr.  (8°) 

where  SVRC  is  the  standard  deviation  of  VRC* 

11.  An  estimate  of  the  kill  probability  of  the  component  by  the 
burst  of  fragments  is  given  by 


PK  =  1  -  e'X  . 


(1  A) 


The  standard  deviation  of  the  kill-probability  estimate  is  given  by 


8PK  =  e"X  8  X  . 


(5) 


Step  11  completes  the  outline  of  a  general  procedure  for 
calculating  a  Monte  Carlo  estimate  of  the  kill  probability  of  a 
critical  component  by  a  burst  of  fragments  where  importance  sampling 
is  used  during  the  picking  of  paths  for  the  sample  fragments.  The 
two  preceding  solutions  are  similar  in  that  sample  fragments  are 
birthed  at  the  burst  origin  and  then  transported  through  any  inter¬ 
vening  inert  material  to  a  possible  impact  with  the  critical  component. 
The  two  solutions  differ  in  that  the  paths  of  the  sample  fragments  are 
picked  from  different  probability  density  functions. 
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c.  The  Adjoint  Transport  of  Fragments  With  Importance  Sampling 


In  the  preceding  method,  we  changed  the  direction  variables  in  the 
kill  integral  (equation  2)  from  6  to  the  coordinates  of  a  point  in 
region  RC  (equation  6E) .  In  the  following  method,  we  will  also  change 
the  remaining  integration  variables  of  equation  6E  from  mQ  and  vq  to 

m^  and  v^.  This  transformation  can  be  regarded  as  interchanging  the 

roles  of  the  source  function  N*  and  the  kill  function  P*.  In  a 
Monte  Carlo  evaluation  of  the  changed  kill  integral,  a  residual  frag¬ 
ment  is  picked  at  the  critical  component  by  sampling  P*  and  then 
transported  back  to  the  burst  origin  (Figure  6) .  The  event  score 
is  then  proportional  to  N*  evaluated  for  the  parameter  values  of  the 
transported  fragment. 


We  assume  that  the  parameters  of  a  source  fragment  are  defined  in 
terms  of  the  parameters  of  the  associated  residual  fragment  by 


mo=  mo  [  n  (r2 )» mi  >v\ ,  8[rz),  g] 


(9A) 


vo  =  vo[n(^)»mi«V|,0(^),‘g] 


(9B) 


r 


2  = 


i 


16 

and  have  continuous  derivatives  over  the  region  of  interest, 
also  assume  that  the  inverse  functions 


(9C) 

We 


=  "1,1/0,  m0,v0 ,0(7^), Is]  , 

(9D) 

il 

_<l 

cF* 

3 

o 

CD 

t£j 

(9E) 

=  ~2  » 

(9F) 

16'W.  Kaplan,  "Advanced  Calculus,"  Addison-Wesley  Publishing 
Company,  Inc.,  1953 
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'  1 


are  defined  and  have  continuous  derivatives  over  the  same  region. 
Equations  9D-9F  are  identified  in  the  integral  of  equation  6E  as 
the  (forward)  transport  operation,  and  equations  9A-9C  will  be 
identified  in  the  following  development  as  the  adjoint  transport 
operation. 

The  kill  integral  now  becomes 


X  = 


r  /•'V’R’tnfe )>^i.V|]Tt[n(^ , g~— To.m0,v0,e(fg)1 
JhcJ°J°  r2  F(r,R) 


N  [  rQ  i  /  vq  <  &  (r2  )]  J  (^i  i  V| )  d  m j  dv| ,  d  V2 


,  (10A) 


J  (  m|(  v, )  = 


d(mo>vo  ) 
d(rnj,  ) 


(10B) 


where  T  [r^  C1^)  >vj ,  g-*-ro ,  mQ  ,VQ  ,0(^3  ]  is  the  adjoint  trans¬ 
port  operation.  The  terms  in  the  integrand  are  rearranged  so  that 
the  source  term  for  residual  fragments  P*  is  leftmost  in  the 
integrand. 

In  its  present  form,  the  kill  integral  in  equation  10A  is  not 
tractable  to  a  Monte  Carlo  solution.  Applying  the  procedure  out¬ 
lined  in  Appendix  A  and  used  in  the  preceding  forward  solution, 
the  integral  is  changed  to  the  more  tractable  form: 


m 


H(r2 )r  r  Umi/V|/ r2'r2') 


RC  RCO  0 


RCOO  r2F(*,R) 

TfR  ( Tz ) (  m,  ,7X ,  9~7q  ,  mQ ,  V0 ,  e(7z )] 


N*Ffr,m0<v0 ,8[rz  )]  J(m|tv,)dm,  dv,  dVg 


1  «  V? 

V  •  ("A) 

VRC 


H(^)=^R*R(^).mMV']drinidv|  * 


(UB) 


L(m,,v|(  r2 ,  rj  )  = 


R*[r,(r2  ),  m)(v,  ]  S(r2  -  r2‘) 


QO  <D 


flp*[^(^)'mi'vi]dnnidvi 


mo 


"cro 


4UL(mi,V,'r2'r2)dm‘dV,dV2=  ]*  niD) 

This  form  of  the  integral  can  be  evaluated  by  applying  the  procedure 
outlined  for  the  preceding  solution.  However,  the  reader  should  note 
that  the  sample  events  used  here  are  strictly  calculational  in  their 
nature  and  do  not  correspond  to  any  actual  ballistic  events.  A  brief 
outline  of  such  a  solution  is  given  below: 


1.  Construct  an  RPP  which  barely  encloses  the  critical  component. 

2.  Pick  a  jaaple-fragment  Fat^  us*n*  method  outlined  in 

the  preceding  forward  solution.  "  1 


3.  Calculate  the  quantities  r^,  IT  ,  and  9^. 

4.  Pick  a  birth  set  of  values  for^the^  mass  and  velocity  of  the 
residual  fragment  by  sampling  L(mj  ,  v^  ,r2»rp  *  The  samPle  fragment 
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is  identified  as  [ (i^)  ^ ,  (m^)  (Vj) .  The  fragment  is  assumed  to  be 

birthed  on  the  outer  surface  of  the  critical  component  which  is  nearest 
to  the  spall-hurst  origin. 

5.  Transport  the  sample  fragment  through  any  intervening  inert 
material  to  the  spall-burst  origin. 


Calculate  the  score  of  the  event  using 


Sr 


ri2  F  ( r; ,  R-  ) 


(12) 


A  fragment  after  transport  whose  mass  (mQ)^  or  velocity  (v  does  not 

fall  within  the  region  of  non-zero  values  of  B*  is  assigned  a  score 

of  0 . 

7.  Repeat  steps  2-6  until  I  events  have  been  conducted. 

8.  Calculate  an  estimate  of  X  using 


_  I 


(VrcISJ/i^VrcS 


(86) 


9.  Using  the  appropriate  forms  of  equations  4C,  calculate  the 
standard  deviation  of  V  _  and  S.  Calculate  the  standard  deviation  of 
X  using 

8X  =  [(\£C8$)Z  *  (SSVff2-  (8D) 


10.  Calculate  an  estimate  of  the  kill  probability  of  the  critical 
component  using 


PK  =  1  -  e'x  , 


(1A) 
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11.  Calculate  the  standard  deviation  of  P  using 


8PK  =  e'XSX  # 


(5) 


Step  11  completes  the  outline  of  a  general  procedure  for 
calculating  a  Monte  Carlo  estimate  of  the  kill  probability  of  a  com¬ 
ponent  by  using  the  adjoint  transport  of  fragments.  An  adjoint 
solution  of  the  test  problem  is  described  in  the  next  section  of 
this  report. 


III.  THE  TEST  PROBLEM  AND  ITS  SOLUTIONS 
A.  The  Test  Problem 

A  cross  section  of  the  test  problem  is  illustrated  in  Figure  7. 

A  spherical  critical  component  having  a  radius  of  10  inches  is  located 
JO  inches  from  a  point  spall  burst.  The  burst  has  an  axis  of  symmetry 
which  lies  along  the  line  connecting  the  burst  origin  to  the  center  of 
the  critical  component.  An  inert  component  of  varying  thickness  (not 
shown  in  Figure  7)  lies  between  the  burst  origin  and  the  critical 
component.  The  inert  component  is  positioned  so  that  it  degrades  or 
defeats  all  fragments  before  they  impact  the  critical  component. 

The  burst  is  described  by  a  truncated  bivariate  normal  distri¬ 
bution  of  fragment  mass  and  velocity.  The  non-zero  values  of  the 
distribution  are  given  by 

B  (r0,m0,v0>  6  )* 


25  N(0)  sin  6  Exp[-0.5Z  0.5Z#] 

(0.9973)*  M(0)  V(0) 


(I3A) 


Z,=  [mo-M(0)]/O.2M(e)  ,  (13B) 


Z2*[vo-V(0)]/O.2  V{0) 


(13C) 


0.4M(S)<  m0  s  1.6MIS)  ,  (13D) 

0.4 V(S)  <  v0  <  1.6V(8)  ,  (13E) 

where  M(0)  is  the  mean  mass  and  V(0)  is  the  mean  velocity  of  fragments 
emitted  at  angle  0  with  the  burst  axis.  The  mean  mass  is  given  by 

M[0 )  =  500-  30000/ it  ,  (T3F) 

and  the  mean  velocity  is  given  by 

V(9)  =  500  +  540003/,  .  (13G) 

The  function  N(0)  is  given  by 

N(0)  »  5  +  (80  9/w  ,  (I3H) 


The  star  superscript  is  dropped  from  the  problem  dependent  quantities 
is  this  example.  We  assume  that  no  correlation  exists  between  the  mass 
and  velocity  of  fragments. 

The  distribution  of  inert  material  between  the  burst  source  and 
the  critical  component  is  described  by  a  truncated  normal  distribution 
at  each  emission  angle  0.  The  probability  P(t)  of  a  thickness  is  given 
by 


P(t )  *  5  Exp[-0.5  Z*,]  /  [o.9973T(0)(2  *)/f],  (131) 
T(0)  *  0. 3  -1.2  9/w  ,  (13  J, 

Z3*  [t-T(0)]  /  [o.2  T(0)] 


(I3K) 


OAT  (6)  s  t  s  1.6  T  (0)  , 


(I3L) 


where  T(8)  is  the  average  thickness  of  inert  material  at  angle  6. 

The  quantity  P(t)  is  used  only  in  the  preceding  equation  to  identify 
the  material  thickness  probability,  and  should  not  be  mistaken  for 
the  kill  function  P*  which  is  used  elsewhere  in  this  report. 

The  kill  function  of  the  critical  component  is  given  by 

R(7J(mp"v, )  =  1CJ8[20m,-*‘ vj  +m|7,][30-D(8)]  (I3M) 


where  D  is  the  distance  from  the  burst  axis  to  the  point  of  impact 
with  the  critical  component  (Figure  7).  The  use  of  the  average 
quantities  and  v^^  would  be  expected  to  introduce  an  error  into 

actual  vulnerability  calculations,  but  we  will  ignore  this  effect 
here  as  it  is  being  investigated  in  a  current  BRL  study?-? 

The  THOR  perforation  equations  (equations  3A  and  3B)  are 
used  to  describe  the  transport  of  fragments.  The  material  coefficients 
used  by  these  equations  are  given  in  statements  30  and  35  of  computer 
program  PR0B3  (Appendix  D)  and  are  not  repeated  here.  All  fragments 
are  assumed  to  be  normally  incident  on  the  inert  component. 


B.  The  Forward  Solution  Without  Importance  Sampling 


The  BASIC  computer  program  PR0B1  (Appendix  B)  is  used  to 
calculate  the  kill  probability  of  the  critical  component  when 
importance  sampling  is  not  used  during  the  picking  of  sample- 
fragment  paths.  The  procedures  developed  and  outlined  in  Section 
II  are  used  here  with  minor  modifications.  Much  of  the  calculational 
procedure  outlined  below  is  given  in  Section  II,  but  is  repeated  here 
for  the  sake  of  completeness. 
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1.  Pick  the  emission  angle  0.  of  the  sample  source  frag¬ 
ment  by  sampling  the  PDF  ®(0)  1 

®{9)  =  2ir{* J* B(m0,v0 ,0)dmo  dv0  =  2ir  sin  0N{0)  .  (14) 
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11  e  picked  values  lie  between  0  and  u/2.  The  rejection  technique 
is  used  to  perform  the  sampling.  The  0-symmetry  of  the  problem  per¬ 
mits  the  calculation  to  be  performed  in  the  xz  -  plane. 

2.  Calculate  the  average  fragment  mass  M.  (equation  131) 
and  t lie  average  fragment  velocity  V.  (equation  13G)  at  emission  angle 

f'l 

i  * 

3.  Pick  a  mass  (mQ).  and  velocity  (v  ).  by  sampling  the 

appropriate  normal  distribution  contained  within  equation  13A.  The 
rejection  technique  is  used  to  conduct  the  sampling. 

4.  Calculate  the  average  thickness  of  material  lying 
between  the  burst  origin  and  the  critical  component  at  angle  6. (equa¬ 
tion  l'Jj.  Tick  a  thickness  t.  by  sampling  P(t)  (equation  131). 

The  rejection  technique  is  used  to  conduct  the  sampling. 

5.  Calculate  the  transport  of  the  sample  fragment  to  a 
possible  impact  with  the  critical  component.  The  residual  fragment 
at  impact  with  the  critical  component  is  identified  as 

o  .  Calculate  the  score  (kill  probability)  of  the  event 
as  P  [  (r  j )  ^ ,  (m^  )  ^ ,  (v^j .  |  (equation  13M).  The  quant  lty  I),  is  derived 

-V 

from  (fj)..  The  score  is  identified  as  S. . 


Y.  A.  Chreider,  "The  Monte  Carlo  Method Pergancn  Press,  Un.. 
island  City ,  New  York,  1966. 
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7.  Repeat  steps  1  through  6  until  I  events  have  been 

conducted. 

8.  Calculate  an  estimate  of  X  using 

.  i  _  - 

X  *  ( N  2  Si )  / 1  *  N  S  .  (4A) 

r 

As  noted  in  equation  ID,  the  quantity  N  is  the  number  of  significant 
fragments  in  the  burst. 

9.  Using  the  appropriate  form  of  equation  4B,  calculate 
an  estimate  of  the  standard  deviation  of  N  and  S.  Calculate  the 
standard  deviation  of  X  using 

8  X  *  [( S  8  N  )*+  (N  8  S*  ]  #  (40 

10.  Calculate  an  estimate  of  the  kill  probability  of  the 
component  by  the  burst  of  fragments  using 


PK  «  l  -  eX  .  (IA) 

11.  Calculate  the  standard  deviation  of  the  kill  probability 
est imate  using 


8PK  =  e  8X .  (5) 

Step  11  completes  an  outline  of  a  Monte  Carlo  forward  calculation 
of  the  kill  probability  of  the  test  component  where  sample- fragment 
paths  are  picked  from  their  natural  distribution.  The  results  are 
tabulated,  along  with  the  results  obtained  using  the  other  two 
methods  of  solution,  in  Table  1  at  the  end  of  Section  III. 
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c.  The  Forward  Solution  With  Importance  Sampling 


The  BASIC  computer  program  PR0B2  (Appendix  C)  is  used  to  cal¬ 
culate  the  kill  probability  of  the  test  component  by  using  the  for¬ 
ward  transport  of  fragments  and  importance  sampling  in  the  picking 
of  sample-fragment  paths.  The  procedures  developed  and  outlined  in 
Section  II  are  used  here  with  minor  modification.  The  calculational 
procedure  used  in  PR0B2  is  outlined  below: 

1.  Construct  an  RPP  which  barely  encloses  the  critical 
component.  In  this  instance,  the  RPP  is  the  cube  whose  faces  lie 
within  the  planes  x-10,  x=-10,  y=10,  y=-10,  z=10,  z=-10. 

2.  Pick  a  point  with  equal  probability  at  any  location 
within  the  upper  right  quarter  of  the  enclosing  RPP.  Because  of  the 
symmetry  of  the  test  problem,  sample- fragment  paths  constructed 
through  these  points  are  representative  of  fragment  paths  for  the 
total  burst. 


3.  Construct  a  ray  from  the  spall-burst  origin  through 
the  sample  point.  The  ray  is  used  as  a  fragment  path  if  it  inter¬ 
sects  the  critical  component.  Continue  picking  sample  points  until 
an  intersection  is  obtained.  A  running  tally  U  is  maintained  of  the 
sample  points  and  is  used  at  the  completion  of  all  sampling  to 
calculate  an  estimate  of  the  volume  VR(.. 


4.  Calculate  the  fragment  emission  angle  9  . 

5.  Calculate  the  average  mass  NT  (equation  13F),  average 

velocity  V.  (equation  13C),  and  average  material  thickness 

(equation  13J)  at  angle  6..  Pick  a  mass  (m  ).,  velocity  (v  ) . ,  and 

1  0  1  0  1 

material  thickness  t .  by  sampling  the  appropriate  normal  distribution. 

6.  Transport  the  sample  fragment  to  a  possible  impact 
with  the  critical  component.  Assign  a  score  of  0  and  go  to  step  9 
if  the  fragment  is  defeated  during  transport. 


7.  Calculate  the  coordinates  of  the  point  of  intersection 

of  the  fragment  path  with  the  RPP  and  with  the  critical  component. 

Calculate  I).  ,  r.  ,  and  R .  . 

i  i  i 

8.  Calculate  the  score  of  the  event  using 


1  • 

[vi]i}G[(r2)i] 

:h 

,R 

i] 

(8A) 
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9.  Repeat  steps  2  through  8  for  a  total  of  I  events. 

10.  Calculate  an  estimate  of  X  using 


X  '(vjbd/i*  VRCS  , 

(8B) 

Vrc  =  (I  V)  /u  . 

(8C) 

11.  Calculate  the  standard  deviation  of  X  using 


—  —  —  £  —  —  2  |/9 

SX=[(VRCSS)  +  ( S  SVRC)  ]  .  (8D) 


12.  Calculate  an  estimate  of  the  kill  probability  of  the 
component  using 


PK  =  1  -  e"X  . 


0  Al 


13.  Calculate  the  standard  deviation  of  the  kill  prob¬ 
ability  estimate  using 


8PK  = 


(5) 


Step  13  completes  an  outline  of  the  calculation  of  the  kill 
probability  of  a  component  by  a  burst  of  fragments  where  fragment 
paths  are  picked  using  importance  sampling.  The  results  are 
tabulated,  along  with  the  results  obtained  using  the  other  two 
methods  of  solution,  in  Table  I  at  the  end  of  Section  III. 
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D.  The  Adjoint  Solution  With  Importance  Sampling 

The  BASIC  Computer  Program  PR0B3  (Appendix  D)  is  used  to  calculate 
the  kill  probability  of  the  component  by  using  the  adjoint  transport 
of  fragments  and  importance  sampling  in  the  picking  of  rays. 

Values  for  the  mass  and  velocity  of  a  sample  residual  fragment 
are  not  picked  by  sampling  P  as  described  in  Section  2C .  Instead, 
the  integral  in  equation  11A  is  changed  to  the  sum  of  K  new  integrals 

where  the  integration  range  of  each  k^  new  integral  is  taken  over 

an  equal  part  of  the  original  integral.  The  total  m^  and  v^  range 

of  integration  is  assumed  to  extend  from  (0,0)  to  the  largest  mQ 

and  vq  values  at  which  N  is  zero.  This  new  representation  of  the 

kill  integral  is  given  by 


k  f  J  dm,  dv,  d  V2 

X  =  if  J  f  - g - 1 - 1 - i 

r  F  ( r  ,  R  ) 


(15) 


where  A.  is  the  integration  range  of  the  m.  and  v  variables  in  the 

tH  ^  * 

k  new  integral.  Equation  15  is  made  less  cumbersome  by  not  including 

the  parameter  lists. 

This  new  formulation  permits  a  set  of  residual  fragments  to  be 
picked  for  each  sample  path.  Each  fragment  in  the  set  is  transported 
back  to  the  burst  origin  and  its  score  is  calculated.  The  scores  of 
all  fragments  in  a  set  are  totaled  to  obtain  the  event  score.  A  step- 
by-step  outline  of  this  calculational  procedure  is  given  below: 

1.  Construct  an  RPP  which  barely  encloses  the  critical 
component.  In  this  instance,  the  RPP  is  the  cube  whose  faces  lie 
within  the  planes  x=10,  x=-10,  y=10,  y=-10,  z=10,  z=-10. 

2.  Pick  a  point  with  equal  probability  at  any  location 
within  the  upper  right  quarter  of  the  enclosing  RPP.  Because  of  the 
symmetry  of  the  test  problem,  sample-fragment  paths  constructed 
through  these  points  are  representative  of  fragment  paths  for  the 
total  burst. 

3.  Construct  a  ray  from  the  spall -burst  origin  through 
the  sample  point.  The  ray  is  used  as  a  fragment  path  if  it  inter¬ 
sects  the  critical  component.  Continue  picking  sample  points  until 
an  intersection  is  obtained.  A  running  tally  U  is  maintained  of  the 
sample  points  and  is  used  at  the  completion  of  all  sampling  to  cal¬ 
culate  an  estimate  of  the  volume  V^. 


40 


I 


picked  ray. 


4.  Calculate  the  fragment  emission  angle  0^  for  the 


5.  Calculate  the  average  mass  NH  (equation  13Fj  and 
average  velocity  Yh  (equation  13G)  of  a  source  fragment  at  angle  6^  . 
Calculate  the  average  inert  material  thickness  t^  (equation  13J)  at 

angle  0^. 

6.  Calculate  the  coordinates  of  the  point  of  intersection 

of  the  fragment  path  with  the  RPP  and  with  the  critical  component. 

Calculate  D. ,  r. .  and  R. . 

i  i  l 

7.  Pick  a  sample  value  for  the  mass  and  velocity  of  the 
t  h 

k  residual  fragment  with  equal  probability  at  any  point  within  the 
integration  range  A^.  This  fragment  is  identified  as 

[(?2)i’(i"l)ik’(7l)ik]' 

8.  Transport  the  fragment  back  to  the  burst  origin  using 
the  adjoint  transport  operation  T^.  As  noted  in  equations  9A  through 
9C,  the  inverse  of  the  forward  transport  equations  would  ordinal ily 
be  used  to  perform  this  operation.  However,  these  inverses  are 
difficult  to  obtain  explicitly  for  the  THOR  equations  3A  and  3B 

and  an  alternate  method  is  used.  First,  we  combine  the  two  THOR  equations 
to  obtain  the  monotonical ly  increasing  function 


9 

=  riv5,)«;"i’cl*scl 

m,=  L  „a ,  TIS  J 

10  (aT) 


C6  CT  CO 

10  (aT)  v0 


2C7+3C8 
2C2+  3C3 


1  r(Vo-v,>Vo  12C2* 
J  L  ,-CI  ,C2  J 

10  (aT) 


(I6A) 


AT  *  aTm0  , 


(I6B) 


A  bar  is  placed  over  the  THOR  material  coefficients  to  clarify  their 
identification  as  single  quantities,  This  use  of  the  bar  should  not  be 
confused  with  its  use  elsewhere  of  specifying  mean  values. 


'j 

r, 


We  then  conduct  a  binary  search  over  the  range  of  possible  values 
to  obtain  a  value  of  (v  )  to  a  specified  accuracy.  The  associated 

(m  )..  is  then  evaluated  by  using 

O  IK  J 


/  -  x  "CS 

(Vq-V,)  Vq 

I0CI  (aT)C2 


3 

2C2+3C3 


(16C) 


9.  Calculate  the  Jacobian  for  the  event.  The  Jacobian 
used  here  is  given  by 


d(m,  ,vt) 

<3(m0,v0) 

and  is  related  to  the  Jacobian  used  in  equation  1 1 A  by 


(I7A) 


jJ  =  l  . 


(I7B) 


The  components  of  the  Jacobian  are  derived  to  be 


dm 

i 

dm0 


_  _  2C7+  3C8 

( 2C7+3C8)10C6(aT)C7vcom  3 

_ o  o 

3m0 


d  m(  -  C6  C7  CO-t  2C7 ♦  3C8 

T“  =  CO  10  (aT)  v„  m,  3  , 

0  vo 

dV,  _  _  Ci  M  M  2C2»3C3 

- =  -{2C2+3C3H0  (aT)  vn  m0  3 

dm0  g 


dv,  —  ci  C2  cs-i  2C2  303 
—  =  I  —  C 5  10  (aT)  v  m  3 
dv0 


(1 8  A) 


( 18  B) 

(18C) 


(18D) 
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10.  Calculate  the  score  for  the  fragment.  The  score  is 

given  by 


Sa*  * 


Am,  Av,  R 

B  ■ 

D^,(m,)iii,(v1kk 

N 

(moLMVoUk,  &»] 

j 

(m0)ik,(v0bk] 

rjF(ii,R;) 

(19) 


where  Am^  and  Av^  are  the  widths  of  the  region  A^. 

11.  Trends  are  identified  during  the  sweep  over  fragments 
associated  with  a  particular  path  which  can  be  used  to  predict  zero 
scores  for  the  remaining  fragments  associated  with  the  path.  If 
non- zero  scores  are  possible  for  the  remaining  fragments,  return  to 
Step  7  and  conduct  the  calculation  for  a  fragment  picked  from  the 
next  Aj..  If  zero  scores  are  predicted,  continue  to  the  next  step. 

12.  Total  the  scores  for  all  fragments  associated  with 
the  path  to  obtain  the  event  score  using 


Si  = 


K 


2S;k  . 


(20) 


13.  Pick  a  total  of  I  paths,  where  each  path  has  an 
associated  set  of  sample  fragments,  and  calculate  the  total  score 
for  each  path. 


14.  Calculate  an  estimate  of  \  using 


x  =  (  ^ )  / 1  s  Vrc  $  * 


(8B) 
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15.  Using  the  appropriate  forms  of  equation  4C,  calculate 
the  standard  deviation  of  VRC  and  S.  Calculate  the  standard  deviation 
of  X  using 


SX^IV^SsfMSSV,*)2]''2.  (8D) 

16.  Calculate  an  estimate  of  the  kill  probability  of  the 
critical  component  using 

PK=l-e'X.  (I  A) 

17.  Calculate  the  standard  deviation  of  P  using 

8PK=  e"XSX  .  (5) 


Step  17  completes  the  outline  of  a  Monte  Carlo  adjoint  calcu¬ 
lation  of  the  kill  probability  of  the  test  component  by  a  burst  of 
fragments.  These  results  are  tabulated,  along  with  the  results 
obtained  using  the  other  two  methods  in  Table  1  in  Section  HIE. 


E.  Test  Problem  Results 


The  BASIC  Computer  Programs  PR0B1,  PR0B2 ,  and  PR0B3  were  installed 
on  a  WANG  2200  computer  and  used  to  calculate  the  kill  probability  of 
the  test  component  by  a  spall  burst.  A  sufficient  number  of  fragment 
histories  were  calculated  to  obtain  a  small  standard  deviation  for 
each  answer.  The  answers  and  their  standard  deviations  are  given 
below : 
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TABLE  I.  Test  Problem  Results 


METHOD 

SURVIVAL  RATE 

NO. 

HISTORIES 

1. 

Forward  Transport  with¬ 

0.3641 

+  0.001 

7000 

out  Importance  Sampling 

2. 

Forward  Transport  with 

0.3644 

+  0.003 

5000 

Importance  Sampling 

3. 

Adjoint  Transport  with 

0.3724 

+  0.014 

2000 

Importance  Sampling 

IV.  CONCLUSION 

The  procedures  for  calculating  the  kill  probability  of  a 
component  by  using  three  different  sampling  methods  have  been  outlined. 
The  defining  equations  have  been  derived  for  each  method.  A  test 
problem  was  devised  and  is  solved  using  each  method.  The  answers 
obtained  by  the  three  methods  agree  within  the  standard  deviations  of 
the  calculations. 

Each  of  the  preceding  forward  Monte  Carlo  methods  is  more  efficient 
for  a  particular  type  of  problem.  The  first  method  is  more  efficient 
for  those  cases  where  fragments  from  the  burst  impact  the  critical 
component  with  a  high  probability.  The  second  method  is  more  efficient 
for  those  cases  where  fragments  from  the  burst  impact  the  critical  com¬ 
ponent  with  a  low  probability.  The  installation  of  both  methods  in  a 
vulnerability  methodology,  where  the  most  efficient  method  is  selected 
for  a  particular  burst  and  component,  should  be  considered.  An 
approximation  of  the  second  method  is  apparently  used  in  the  present 

version  of  the  VAST  methodology . 7’ ^  The  third  method  (adjoint 

transport  of  fragments  with  importance  sampling)  is  not,  to  the 
author's  knowledge,  used  in  any  vehicle  vulnerability  methodology. 

In  many  vulnerability  studies,  a  comparison  of  the  kill  probability 
of  a  particular  vehicle  by  the  spall  from  different  munitions  is  desired. 
'Phis  comparison  is  presently  obtained  by  calculating  the  disabling  prob¬ 
abilities  of  spall  bursts  produced  by  the  different  munitions  for  a 


D.F.  Haskel  and  M.J.  Reisinger ,  "Armored  Vehicle  Vulnerability  Analysis 
Model  -  First  Version,  Introduction,"  USARRADCOM/BRL  Report  No.  1857, 
February  1976. 
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common  set  of  spall-burst  origins  and  a  common  set  of  sample-fragment 
rays  from  each  origin.  A  detailed  description  of  the  passage  of  each 
ray  through  the  vehicle  is  calculated  and  stored  on  magnetic  tape. 

The  loss  of  capability  in  a  component  due  to  a  particular  munition  can 
then  be  obtained  in  a  second  calculation  by  using  the  retreived  ray 
descriptions  in  conjunction  with  the  spall-burst  source  term  of  the 
munition. 


A  comparison  of  the  kill  probability  of  the  spall  bursts  from 
different  munitions  can  also  be  calculated  using  the  adjoint  transport 
of  fragments.  An  adequate  number  of  sample  Tesidual  fragments  are 
picked  and  transported  to  define  the  adjoint  fragment  flux  at  the 
burst  origin.  This  adjoint  flux  can  then  be  convoluted  with  any  spall - 
burst  source  term  to  obtain  its  kill  probability.  An  advantage  of 
the  adjoint  method  over  the  forward  method  is  that  less  information 
is  stored  for  the  second  calculation  and  the  second  calculation  is 
shorter.  Additionally,  in  the  adjoint  method,  sample  fragments  can 
be  picked  more  efficiently  to  obtain  a  specified  statistical  confidence 
level  for  the  calculated  kill  probabilities  of  a  large  number  of 
different  burst  descriptions. 
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APPENDIX  A.  THE  MONTE  CARLO  EVALUATION  OF  DEFINITE  INTEGRALS 


The  Monte  Carlo  method  can  often  be  effectively  used  to  evaluate 
complex,  multi-dimensional  definite  integrals  which  aren't  amenable  to 
evaluation  by  other  methods.  This  method  is  outlined  below  for  a  two- 
dimensional  integral,  and  the  terminology  commonly  used  is  defined. 

The  outlined  method  can  be  easily  extended  to  integrals  of  a  higher 
order  of  dimensionality  by  iterating  the  outlined  procedures. 

We  wish  to  evaluate  the  integral 


Xs/*  f  F (x  ,y)  Q(x,y)dydx  (A-|) 

xi  yi 


where  F(x,y)  and  Q(x,y)  are  continuous  and  single  valued  within  tne 
specified  region  R  of  integration.  This  integral  is  normally  inter¬ 
preted  as  the  volume  under  the  function  F(x,y)Q(x,y)  which  is  located 
within  the  integration  limits.  However,  if  we  normalize  Q(x,y)  over 
R  to  1,  equation  A-l  becomes: 

xi  y* 

X*C  f  f  F(x,y)  P(x,y)dydx  ,  (A-2) 

x«  yi 


Vji 

C*/  */  *Q(x,y)dydx  ,  (A-3) 

xi  y, 


P(x,y)  3  Q(x,y)  /  C  ,  (A-4) 


and  the  integral  in  equation  A-2  can  now  be  interpreted  as  the  average 
expected  value  of  F(x,y)  over  R  where  the  values  of  x  and  y  are 
predicted  by  the  probability  density  function  (PDF)  P(x,y). 


47 


The  quantity  X  can  now  be  easily  estimated  by  using  the  second 
interpretation  of  the  integral.  Several  sets  of  values  are  picked 
for  x  and  y  by  sampling  P(x,y),  and  the  mean  value  of  F(x,y)  is 
calculated  for  these  values.  An  estimate  of  X  is  given  by 

Xs  C  2  F{ y^)  /I  *  C  F  (A- 5) 


where  i  is  the  sample  index  and  I  is  the  number  of  sample  sets  picked. 

A  bar  is  placed  over  a  quantity  to  identify  its  mean  value.  The 
quantity  F(x,y)  is  often  identified  as  the_estimator ,  and  F(x-,yi) 
is  identified  as  the  score.  The  quantity  X  converges  toward  £he  true 
value  of  X  as  the  number  of  sample  events  is  increased. 

19 

Various  sampling  procedures  are  available  for  picking  sample 
values  for  a  variable  (variant)  from  a  one- di mens ional  PDF.  Therefore, 
we  will  change  equation  A-2  to  a  form  where  the  two-dimensional  PDF 
P(x,y)  is  replaced  by  a  pair  of  one-dimensional  PDF's.  This  change 
is  given  by: 

x'  xx 

X  =  c  /  *{/ jf  FU.yJYtx.x'.yJdydx}  X(x')dx'  ,  (A-6) 

X1  x,  y. 


.y* 


X(x')  =  /  P(x',y)dy  , 


y. 


(A-7) 


P(x,y)8(x-x') 

y* 

f  P(x',y)dy 

y, 


(A- 8) 


where  6(x-x')  is  the  Dirac  delta  function  and  x'  is  a  dummy  value  of  x. 


E.D.  Cashwell  and  C.J,  Everette ,  ", A  Practical  Manual  on  the  Monte  Carlo 
Method  for  Random  Walk  Problems Pergamon  Press  (1959). 
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We  will  now  outline  a  step-by-step  procedure  for  calculating  a 
Monte  Carlo  estimate  of  X.  This  procedure  is: 

1.  Pick  a  value  of  x'  by  sampling  X(x').  This  value  is  identi¬ 
fied  as  xf . 

1 

s 

2.  Pick  an  associated  value  of  y  by  sampling  the  one-dimensional 
PDF  Y(x,x',y).  This  value  is  identified  as  y' . 

3.  Calculate  the  score  F(x'^,yj)  for  the  event.  This  quantity 
is  identified  as  F.. 

4.  Calculate  the  scores  for  a  total  of  I  events  by  repeating 
steps  1-3. 

5.  An  estimate  of  X  is  calculated  as 


X  =  C  2  Fx  =  C  F  . 


(A- 5) 


6.  A  measure  of_the  confidence  level  of  X  is  given  by  the 
standard  deviation  of  X.  This  quantity  is  given  by 


X  2  —  *  1/ 

r  X(CF*)-IX  1  '* 

X  =  L  I(I-l)  J  * 


( A-9) 


Step  6  completes  a  step-by-step  outline  of  the  Monte  Carlo 
evaluation  of  the  integral  in  equation  A-2.  Often  an  integral 
is  more  amenable  to  evaluation,  or  the  statistical  convergence  of 
X  can  be  attained  more  readily,  if  the  integrand  is  rearranged  from 


t 
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the  form  in  which  it  normally  occurs.  This  importance  sampling  is 
developed  below  for  equation  A-l: 


x*  y* 


/  /  F(x,y)Q(x,y)dydx  * 
y  i 

;  y  [F<».,)G<x.y>]  dydx. 

x,  y(  ' 

( A-IO) 

x2  ya 

/  /  FU x , y )  q'(  x,y  )dydx  * 
x,  y, 

( A- 1 1 ) 

x 

c'f  7  V(x,y)P'(x,y)  dydx , 
xi  yi 

(A- 12) 

xx  yt 

c'=  f  f  Q'(x,y)dxdy, 
xi  yt 

(A-13) 

p'(x,y)  MQ'(x.y)  /C'. 

(A-14) 

liquation  A- 12  can  now  be  evaluated  by  using  the  same  procedure  used  to 
evaluate  equation  A-2.  However,  the  sets  of  variable  values  (x^y^) 
and  the  scores  F' .  as  well  as  the  normalization  constant  C' 
will  differ  between  the  two  calculations. 


APPENDIX  B.  PROB  1  LISTING 


RtM  F'F’OGR'RM  PPUB1 

10  F  EM  FORliRPO  SPhLL  TEST  PROBLEM 

11  PEM  NO  I  MPOF'TRNC  E  SRMPLING 
1  —  ^ElEC  T  P 

i:  SELECT  PPINT  .ill; 

15  INPUT  "RANDOM  NUMBER  SEED"  ■  ZB 
10  INPUT  " NUMBER  OF  HISTORIES". II 
15  P‘1—0  P2=0  J=0 

10  INPUT  "THOR  VELOCITY  COEFFICIENTS" , Cl, C2,  Cl. C4, C5 

15  INPUT  "THOR  MRSS  COEFFICIENTS". C6, C7, C8, C9, CO 

16  INPUT  "ORTE" ■ XI 
print  "DRTE="  ::i 

IS  PRINT  "NUMBER  OF  HISTORIES*", 11 

40  F1*="C1=" :  F2f="C2=“  F2**"C2=" :  F4*="C4="  F5*="C5=" 

41  PRINT  "THOR  VELOCITY  COEFFICIENTS" 

42  PR  I  NT  US  I  NO  48,  FI*.  Cl,  F2f,  C2,  FI*.  C3,  F4*,  C4,  F5*.  C5 

41  Flf="C6="  F2*="C7="  F3*="C8="  F4*="C9="  F5f="C0=" 

44  PRINTUSING  48.  F‘l*.  C6,  F2*,  C7,  FI*.  C8,  F4*.  C9-  F5*,  CO 

481###  -##  ###  ###  -##.  ###  ###  -##  #4#  ###  -##  ###  #4#  -##  ### 

49  PRINT 

50  B1=I000.  ’I  14159 
55  B 2 =54000/' 1.  14159 
60  81=180/1.  14159 
65  B6=RRCSINOH.  5> 

1 0  B4=  •  £:  1’  ♦  B6+5  >  +  S I N 1  E'6  y 
75  B5=SQP<2+2.  14159 > 

80  B7 =1.  2/2  14159 

99  SELECT  PRINT  005 

100  FOR  I =1T0  II 
105  REM  PICK  R1 
110  r=j+ i 

115  GO  SUB  “'99 
120  R1=Z9+B€ 

125  GOSUB  ' 99 

120  X1*Z9*B4 

125  X2= < B3*Rl+5 > *S I N < R1 > 

140  IF  X2CX1THEN  110 
145  PEM  PICK  M 
150  M O = -  B 1  *  A 1 + 5 0 5 
155  GOSUB  1 0  •-  M 0  > 

160  M=X2 

165  REM  PICK  V 

1 7  0  V  0 = B  2  *  A 1 + 5  0  0 
175  GOSUB  1 0  < V O  > 

180  V=X2 

180  v=::2 

185  REM  PICK  T 
190  T 0=-B7+Rl+0.  2 
195  GOSUB  10 ( T 0  > 

200  T  =X2 

205  PEM  CRLCULRTE  TRANSPORT 


SI 


.'C  GOi-UB  TO 
.10  GOSlJB  .201  1  1 
.15  IF  M  >3 THEN  .227 
12 5  OOS'JB  20'  2  • 

1 ... 5  IF  OTHEN  270 
_  p~ m  GOTO  25*3 
ii  ijijf.ijB  EE' 

.  _'S  r  -  •  20  ♦  M+V-t-M  +  V  1  +:  1  EO-P  1  '1+DO00OQ00 
.5  I-  P  11THEN  245 

..40  P-1  0 

.45  P1=F'1+P  P2=F'2+F'+'P 

.5m  PRINT  I  ■  T.  P 
1 " 5  NEXT  I 

ELECT  PRINT  215 
'•<  P=Pl.  Il 

"5  P2=-  P2-I1+P+P  11+  '  1 1-1)  • 

=  2+:#F'I  +: I  1 +  B4  +  EE.  T 

2= ■  [i+  j-ii+u  j  ♦  ;+•  T-i  ■  •. 

2=22+' \ 

.."r+O  P2=*  F'+:22  .« "'2+'  P2  +  :  ;  •  "2 
725  P2=S0P 1  P2 > 

. 00  F'-F'  +  2  P— Ei'F''.  — F  ) 

.  05  =,L=F"+:P2 

1.0  PRINT  "SUPVIVIHL  F'FiTE  =".P 
111  PRINT  "S  P  STRNORPO  DEV  I RT  I  ON  =".P2 
15  Ef 4D 

■ 'ML1  DEFFN  20 1  7.1' 

_m05  F  EM  CRLCULRTES  TF'RNSF'OPT'  THOP  ' 

1010  IF  7171THEN  2025 
1014  Ml=f'1 

15  M*M— •  10  *  Co  >  + '  T  +'R0  "C7  '  *•  M'“'C8  +m  V'CCO 
.020  GOTO  2020 

1025  —  •  10* Cl  '•■+•'  1  T  +  RO  •  "C2  '  +  1  M1*  C2  '•  +  L  VC5 

IMIO  RETURN 

.poo  DEFFN  20 

1.505  REM  CRLCULRTES  RO 

.5  10  L=M.  '14  42  2 


5110  ;  !4=5  0*<'X2-X1>/X1 
5125  X4=EXP*-0  5*X4~2.> 

5140  X4s5*X4/'<  0.  972+B5+X11 

5145  Tc  '4  XX 2 THEN  5110 

5150  RETURN 

SiOt1  DEFFN  68 

6205  REM  C0LCUL0TE  R 

6210  X1»-T0N<.01> 

6215  09-1+X1+X1 
6220  B  9 = -  2 + X  2 + X 1  +•  X 1 
6225  09=  X1*X2 >  '~2— 100 
6220  R=69*69— 4*09*09 
6225  IF  R20THEN  6250 
6240  PRINT  "NEG0TIVE  R",R 
6245  GOTO  6265 
6  ^  5  ‘J  R  “  ^\r.tp  v  R'  1 
6 2 5 5  R=0  5* 1  P- B 9 > / H 9 
2  K'  u  R — > ,  i  +-  *  p— 2  y 
6265  RETURN 
7000  DEFFN ■ 77 
7005  REM  DEBUGGING  SR 

7010  H1*="I=":  02*="J=“ •  02*= "01*" :  04*= "R=" •  05*="P1=" 
7015  06*="R2=":  07*="M="  :  08*="V="  :  89*=’^="  00*="?'=" 

7020  PRINTUSING  7100,  Ul*,  I,  02*.  J,  02*,  01,  04*,  R,  05*,  R"1 
7025  PRINTUSING  7100,  06*,  R2,  07*,  M,  08*,  V-  09*,  T,  00*.  P 
7100X444  4444  44  444  4444.  44  444  4444.  44  444  4444  44 

4444  44 
7105  RETURN 
9W0W  DEFFN  99 

9005  REM  P0NDOM  NUMBER  M0IN  SR 

961101  Z8=25*Z8 

9015  GOSUB  98 

9028  Z8=5*Z8 

96'25  GOSUB  ''  '98 

9 0 2 0  Z  9 = Z  8  .-'671 0 8 86 4 

9025  RETURN 

9050  DEFFN ' 98 

9055  REM  P0NDOM  NUMBER  SECOND0PV  SR 
9 0 6  0i  X  5 = Z  8  /  6  7 1 0i 8 8 6 4 
9065  X6=67 1018864+  I  NT  <  K5  > 

96'76’  Z8=Z8— X6 
9075  RETURN 
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APPENDIX  C.  PROB  2  LISTING 

1  POT  P'PiHiPAM  PRl.lB  X 
.  I  ■  I l  l  SI  S2<'2'»..  W<3> 

Jh  REM  FORWARD  SPALL  TEST  PROBLEM 

It  REM  NO  IMPORTANCE  SAMPLING 

!  SELEC  T  P 

1  .  SELEC  r  PRINT  215 

J.".>  INROT  "RANDOM  NUMBER  SEED".  ZB 

,ii  INPUT  "NUMBER  OF  HISTORIES"..  II 

-"-.i  Pl~0  P2=0  J=0 

'  i  1  NRU  T  "  THOR  VELOC I  TV  COEFF I C I  ENT'S  "  .  Cl .  C2,  C3,  C4 ..  C5 
. INPUT  "THOR  MASS  COEFFICIENTS",  C6- C7,  08 •  C9,  CO 
.  INPUT  "DATE".  XI 
PRINT  "DATE=  ' .  XI 

■  PRINT  "  Nl  JMBER  OF  H  J.  STOP  I  ES=  "  .  1 1 
■4,.i  Fil="r,l*"  •  F2T="C2~"  :  F2f.=  "C2  ="  :  F4*="C4="  :  F5*="C5~" 

•U  PRINT  "THOR  VELOC I TV  COEFFICIENTS" 

4,  PR  IN  riJSING  48-  F:l*,  Cl.  F2*.  C2,  F3*,  C3,  F4f .  C4.  F5f .  C5 

r  i  f  - "i  8="  F2t  =  "C?="  F2*="C8="  F4*="C9="  F5f="C0=" 

44  PRINTUSING  48.  FI*.  C6,  F2t,  C7,  F3t,  C8,  F4*,  C9.  F5$,  CO 
4  -3.444  -44  444  444  -44  444  444  -44  44#  444  -44  444  444  -44  444 
•  I *-'P’  !.  N I 

■':>  u  i-.,  i.  -- 3  0 0 O  X  2  1 4 1 5 9 
r:o  B 2 -7,4000X3  14159 
r-,-:i  8  2~  180X3:  14159 
,■  5  Ei5.i~S.i4R  ■:  2*  ■■■  1415,9  > 

8U  BX=1  2X3  14159 

85  C -10  SI  •  1  1  ~20  '  S  1  •:  2  >  =0  :  S:1  <  3  >  =0 
99  SELECT  PRINT  005 

1.00  For  i— itu  ll 

1 0  ",  , j i. )S  *  I B  '  P. 7 
1.0 7  REM  RIO  M 
1  1.0  M  0  -  -  A 1  ■*  B 1 + 5  0  5 
115  GiV- 1  ip:  10-  j'TO  • 

180  M=0X 
185  REM  RIO:  V 
.l  _,0  V0  =A  1  402+500 
135  ,  ,p  '  10 <  V0 > 

I 4U  V-U2 
1.4  5  RR--1  PI i*l  T 
L5-0  T 0s=— A 1  +.B7 -*-0  3 
185,  i  ] OS  US  lO'  TOi 

180  i=02 

.1.85  REM  CALCULATE  TRANSPORT 
1 t 0  GO SUB  30 
1  •• ' *i>  GO  SUB  20  1  '• 

.1. -ii  IF  MC0THEN  195 

i . .-  >  1 1 1  i  — .  i  I  301  2  ■' 

150  IP  V-OTHEN  200 
1.  95  P'--0  GOTO  235 

200  P  -  <:  20  +  M+V+M+V  +  •:  30-R  >  / 100000000 
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PRESIDING  PAfll  BUWK-WOT  Tll^ED 


n‘.  If-  P  I  THEN  2 IS 
:l‘i  P-1 

:!  I  S  P  -'P  +  '.  B3+H 1+S  • 
t.'l  I-' 1.  ~  '  R2  +  R‘l  1  p1  ,_l 
;’2S  P  :*P,  '‘  1  +  121+01+"  I  • 

_  +1  1'  .  +  F',-'  *:  R  1  +  R'V  +  F',‘  < 

•  S  PRINT  I.  T .  p 
:4o  PL=P1+P  F'2=P2+P+F' 
:4  LSELEL  T  PRINT  CIS 
'4.,:'  fiOSIJB  ■"?? 

;4S  NEXT  I 


Sm 

E‘ 

=Pi/I 1 

-;iS 

p- 

2-' 'P.2-1 1 ,+P +P 

11+ ■ I 1—1 > 

r,  0 

P 

2 -SDR 1  P2 ■ 

r-.  S 

0 

■  11  +  .7-11*  I  IT/ 

•  T+  T  + ■  T-  1 

••0 

o 

2=S0R  •  02 ■ 

,’S 

0 

=  '.’.2+0  1  +  1  L.-'.T 

HI-1 

r.i 

2=02+ 0 

8S 

p 

2- '  P  +  02  1  ■" 2+ 1  O+P 

2  >  ~ 

90 

p 

2= SOP  •:  P2 

;9S  R=P+0 
IS?  P--EXP--  -P"> 
"4:-;  R.>=P+P2 


00  SELECT  PRINT  1 IS 

10  PRINT  "SURVIVIRL  RHTE  •  P 

■  11  PRINT  "S  R  STHNDhRD  DEV  I RT ION  -".PC 

•  IS  END 


TOl'S 


1000 

L'PPPN 

20'  X 1  ■ 

REM  C 

ftLCOLFi  !E 

:0  ip 

IF  2.1 

Tl THEN  2 

014 

Hl-M 

VI  IS 

M=fT-' 

10  CP  ■  + • 

>•1,4-1 

i  TOT  Ti 

2020 

M2  *5 

V=V  -  > 

10  *  C 1  :■  +  ■: 

;020 

RE  TOR 

N 

jShm 

r.'EFPN 

'  20 

REM  C 

HL.C01.hTE 

010 

01=H/ 

14  432 

SIS 

02=01 

/7  764:1 

c.  I  CL  0 

Ij's-i  0 

7S+02/3 

~i  l 

141S'-'  +  02 

H0-HO 

2  S4  +  2 

S  20 

RE  TOR 

N 

loo 

DEF  P  N 

10  •:  01  • 

100 

REM  S 

ftMP'LE  FP 

110 

000 1  IE: 

•-4  C| 

1  IS 

02=01 

+  \0  4+1 

120 

oosue 

9  9 

12S 

•  ■ :  -P:S 

+ol +0  97 

'  C  S 


i+Zs 


12+  i 1  . 


i +  7’+,  I  - 1  - 
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5i20 

04=5.  02— Ql > .'’01 

51  "5 

0  4=EXP  <  -  01  5  +  Q  4 2  > 

5140 

r  4=5  +■  0 4 '  ■  0  973  +  B 5 + 6!  1  > 

5145 

IF  04 70 3 THEN  5110 

5150 

RETURN 

h.MMM 

OEFFN  "67 

6005 

REM  PICK  Ft  NO  DEFINE  S 

HOT  LINES 

if.  007 

J=  T  +  l 

6010 

LjUbUB  ’38 

6012 

S2 1  1  =C*  C  2+Z9-1  > 

6015 

GO  SUE:  •"  99 

6020 

82  v  2  >  =C*:Z3 

6025 

GO SUB  "99 

6030 

S2«  3  =C+:Z3 

6035 

R9  = 1  Sic  1  > -S2 C 1 .)  >  "2+  C S 

1 C 2>— S2 C  2  > > 

6040 

P3=bUF:  •  Ry  •' 

6045 

W '•  1  ■  =  C  S2  <  1  >  -  SI  C 1  .>  > /R9 

6050 

W  '■  2  '■  =  •:  S2  <  2  > —SI  C  2  >  >  /R'9 

6055 

W <  2  >  =  C  S2  C  3  >  —SI  C  ±  >  >  ."’R9 

6057 

GO'S  UE:  "68 

6058 

IF  Fl=— 1THEN  6007 

6060 

xi=C 

6*06*5 

VI =S 1 < 2  > +W < 2 > * < C— SI C 1 

>  ' .WC 1 

6070 

21=S1C  3  > +W  <  3 )  *■  <  C-S1C 1 

>  )/W(l) 

6075 

■-•*  —  ~  1  ' 

6080 

V2=S1  <  2  > -W  C  2  >  *■  c.  C+S1C 1 

))/l.Kl) 

6085 

Z2=S1  C  3  >  — W<  3  >  +■'  C  C+Slc  1 

/W  C 1  > 

6080 

V2  =C 

6085 

>'2  =S1  C 1  >  +W  <  1  >  *  c  C— SI  <  2 

:>  .'W  c  2  > 

6iO0 

Z2=S1 ■ 2  > +W C 3  > *  < C-S1C 2 

>  >/W<2> 

6105 

Z4=C 

6110 

X 4 = S 1 < 1 >  + WC 1 > * < C-S1C 3 

>  )/W(3) 

6115 

V4=S1  2  >  +W < 2  >  *  C C— SI  <  2 

>/WC3> 

6117 

REM  CFtLCULflTES  R2  RND 

Ftl 

6120 

R2=c  SI  C 1  > -XI  ‘2+  <  SI  <  2 

>-Vl>’'2+CSl 

6125 

R2=S0Ri  R2 

6120 

Ftl=FtF'CCOS  <  10.-'R2  > 

6125 

PEN  CRLCULFfTE  RCL 

6  J_  4  ^*.1 

01=  C  si '  1  '  — X2>  '  2+  c  51  2 

>-V2>  "2-KS1 

6145 

02='  SIC  1 — X2 > -2+ SI 1  2 

>  -V3  >  '*'2+  <  SI 

6150 

03=  ■  SI 1  1  >— X4  ’<  "2+  < SI c  2 

>  — V4  >  "'2+  C  S'l 

6155 

04=01 

61601 

IF  021:  04 THEN  6170 

6165 

04=02 

6170 

IF  02 ! 04THEN  6180 

6175 

04=0(2 

6180 

R‘i=bUR '  i2-!,4  ■' 

6185 

R1=R1-R2 

6  180 

RETURN 

6300 

OEFFN  68 
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F205  PEM  CRLCULRTE  P 
F207  F  1=1 

F  2  1 0  A  9 =1+  W  2  1 " 2 + M  <  3  > "  2  >  /  <  W  <  1  >  ■'  2  > 

o2 13  B9=2+W  ■-  2  *  <  SIX  2  > -W <  2  >  +S1  <  1  > /W <  1  >  >  /W  1  > 

F220  B9=B9+2+W-  3  :■•+  0 SI 2 > -W <  3  >  *S1 < 1 > /W  < 1 >  >/W<  1  • 

-223  C9=S1<  2  2— S 1  <  1  1  +  U ■'  2  >  *  (  2 +S  1  <  2  >  —  U C  2  > + S 1 1 . 1  1 X W •  1  1  v'W<  1) 

F 2 1 0  C9=09+S1  •  2  :•  ~2-Sl <  1 >  *W  <  3  >  *  <  2*S1  <  3  > -W  0  3  >  *S1  ••  1>/W<  1  '  >/W*  1  ' 

F2 I 5  09*09—0*0 

-  240  U1=E"9+B9— 4+R9+C9 

F245  IF  0 1>=0THEN  F2F0 

F250  F1=-:L 

F 2 5 5  GOT 0  F  3  4 0 


*260 
*  2l  6  5 

01=SQR<Q1 > 

X5=0i  5*  <.  -B9+Q1  > 

XR9 

t£.  270 

XF=0.  5+  <  -B9-01  ;■ 

XA9 

f_-  c 

V5=S1 0  2  >  +W  0  2  :•  *  *• 

xs-si 1 1  :• 

,‘W 1  1 

*230 

V  F = S 1  •'  2 .)  +  U  •  2  + ' 

KF-Sl 

.'XN'l 

*2*5 

Z5*S1<3 >  +W < 3  >  +  < 

X5-S1 < 1 > 

'/WCl 

*2'?0 

ZF=S1<  3  >+W<  !>*•:. 

XF-S1-:  !  .:• 

■XL'Kl 

F293  Rl=  i:.  SI  ■:  1  >  — X5  •'  "‘2+  < SI <  2  > -V5  > "2+  < SI < 3  > -Z5  >  ‘  2 
—  _ R4  =  < si 1. 1  —  Xo ’  2*1-. '-.I1. 2.,—  VF.-*  S'l1'  2'.'1  — IF  1  2 

FI05  IF  PI 2P4THEN  6315 
FI  10  X5=XF  :  V5=VF  :  I5=ZF 
FI 13  r=v5*V5+Z5*Z5 
F_  20  R'  = t. 1 F' 1  R .) 

FI40  RETURN 

TOGO  OEFFN  "  77 

7003  REM  C'EBUGG  I NG  SR' 

7010  Rif  =  “  I  =  ”  R2-f="  J="  :  RIf="Rl="  :  A4f="R="  R5f="Pi=" 

7015  AFf="R2="  R7f="M="  •  R8f=,"v,=  "  :  fl9**"T="  :  H0f="P=" 

7020  PR  I  NT  USING  7100.  IJlf  •  I,  R2f,  J,  Rif,  Rl.  R4f,  R.  R5f  .  PI 
7025  PRINTUSING  7100,  RFf,  R2,  R7f,  M,  R8f,  V.  R9f,  T.  ROf .  P 
7100V###  ####.  ##  ###  ####,  ##  ###  ####  ##  ###  ####  ##  »#* 

####  ## 

7105  RETURN 
9  0  UU  D  E  F  F  N  ■'  9  9 

9005  REM  RANDOM  NUMBER  MR IN  SR 

90110  18= 25 +19 

90'15  G0SU6  98 

9  0i  2  0i  18=5  +  Z  8 

9023  GOSUB  "'98 

ROC 0i  Z9=Z8.  'F71088F4 

9015  RETURN 

9030  DEFFN'98 

9055  REM  RANDOM  NUMBER’  SEOONDRRV  SR 
9  0i F 0i  05=I8XF7 1 0i 88F4 
90'F5  OF=F71083F4* 1  NT  <05  > 

90i 70  Z8=Z 8-0IO 
-075  RETURN 
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APPENDIX  D.  PROB  3  LISTING 


1  REM  PROGRAM  PR063 

2  DIM  Sit  2  • .  S2v2 >.•  Wv2.:- 

10  PEM  ADJOINT  SPALL  TEST  PROBLEM 

11  REM  NO  IMPORTANCE  SAMPLING 
11  SELECT  P 

I!  SELECT  PRINT  215 
15  INPUT  "RANDOM  NUMBER  SEED",  Z8 
20  INPUT  "NUMBER  OF  HISTORIES", II 
25  P1=G  P2=0:  T=0 

20  Cl=6.  399 :  C2=0.  889 :  C3=-0.  945 :  C4-1  262 :  C5=0.  019 

25  C6=-2.  507 :  C7=0  128 :  CS=0  825 :  C9=0.  142 :  C0=O  761 

26  INPUT  "DATE",  XI 

27  PRINT  "DATE*",  XI 

28  PRINT  "NUMBER  OF  HISTORIES*",  II 

40  Flf="Cl*“ :  F2f="C2="  F2f="C2=" :  F4$="C4=" :  F5f="C5=" 

41  PRINT  “THOR  VELOCITV  COEFFICIENTS" 

42  PRINTUSING  48.  Flf.  Cl,  F2$,  C2,  F2#,  C2,  F4$,  C4,  F5f  •  C5 

42  Flf ="C6=" :  F2f=“C7=" :  F2f="C8=" :  F4f="C9=" •  F5*="C0=" 

44  PPINTUSING  48,  Fit.  C6,  F2t,  C7,  F2f,  C8,  F4f .  C9,  F5f .  CO 

48;:###  -##  ###  ###  -##.  ###  ###  -##.  ###  ###  -##  ###  ###  -##,  ### 

49  PRINT 

50  Ea  =  1000.  '2  14159 
55  62=54000.  '2-  14159 
60  62=188/2  14159 

65  C  1=10 ‘Cl:  C6=10-'C6  :  J0=1 
67  E2=- 2*C?+2*C8>/2 
69  El=- 2*C2+2*C2>/2 

71  E2=E2.  FI 

72  E4=1/E1 

75  B5=S0R  2*2  14159  > 

80  67=1.  2/2  14159 

65  C  =10  .  SI '/l >*20  SI- 2  ,=0  Sl<3>*0 

90  A=4*#F'  I  *14.  422*7.  7641 
92  A=  (' 3/A  666667 

94  A=#PI+A/<2  54*2.  54 > 

99  SELECT  PRINT  005 

100  FOR  1=1 TO  II 

105  REM  PICK  ADJOINT  ORIGIN 

HO  GOSUB  '67 

114  M0=-B1+ Al+500 

1 16  70=62* Al+500 

120  M2=l  6* MO :  M1=0.  4*M© 

140  72*1.  6*70 :  71=0.  4*70 
155  REM  PICK  T 
160  T0=— A1+B7+0  2 
165  GOSUB  10  -.  TO  > 

170  T *02 

175  M9=75  79=1250 

ISO  E5=- A*T  -  - C2 
185  E6=- A*T  >  *‘C-7 
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200  GO  SUB  00 
200  NEXT  I 
210  P=P1/I1 


v  20 

P2= 

•  P2— 

ii 

+F'+F‘  >  /  <.  11+  <  1 1-1  > 

770 

F2= 

!z»WP  •*. 

P2 

140 

Q2= 

<  I  1* 

J- 

11+ 1 1.:-/'  T+T+CJ-l 

1-50 

Q2- 

SOR< 

02 

> 

1 60 

Q=i: 

2*C 

3>+Il/J 

270 

1  >!  ci = 

Q2:*Q 

780 

P2= 

P*Q 

2> 

~2+  <  Q+P2 ')  "'2 

230  P2=SGR < P2  > 

400  P=P+0 
405  P=EXP  — P  :< 

4 OS  P2-P2+P 

410  SELECT  PRINT  215 

420  PRINT  " SURV I V I ML  ROTE  =“,  P 

410  PRINT  "SR  STANDARD  DEVIATION  =  ".P2 

450  END 

1000  DEFFN'OO 

1004  SELECT  PRINT  005 

1005  REM  DIRECT  MASS  AND  VELOCITV  SAMPLING  AND  SCORING 

1006  Q1=C R1+R2 > /R2 

1007  01=R1+  <  1+Q1+Q1+Q1 ) /2 

1008  P9=25+M9+V9+ < B2+A1+5  > / < 2+#P  I  +M0+V0+R2  +P2+0.1  > 

1010  K  =— 1 :  V4=V2:  P€=0 :  F9-100 

1015  Q0SU8  '39  :  Z7=29 
lOlS  GObUB  ''  39  :  2  6  ~  Z  3 
1020  k>K+l :  L=-l 

1025  VS*  <  1+27  >  +V9 

1026  IF  V33V2THEN  1215 
1020  V3=V8 

1040  L=L+1 
1045  MS*  <'  L+Z6  >  +M9 
1 0 5 0  G  Ci  S  U  B  7  7 
1055  V=Q7 

1060  REM  CALCULATE  M 
1 0 6 5  Q 1 = C 1+E5 + < V ' C 5  > 

1 0 7 0  Q  1*<  V — V  8 > /  Q 1 

1075  M = Q 1  >'■.  2  ,•••'  <  2  *  C  2+2+ C  2  >  "> 

1112  01=' R1+R2  > /R2 
1114  Q 1 = P 1 +  ■:  1 +01+ Q 1  +  0 1  >  /  2 
1150  IF  M2M2THEN  1190 
1155  IF  M1M1THEN  1190 
1160  IF  V7V2THEN  1190 

1165  IF  VI VI THEN  1190 

1 166  01=5+  V- VO  >  /VO 

1167  01=0  5+01+01 
1163  F‘S=EXP •!  -01  > 

1 1  ^0  01=5+  v  M-MO  /MO 
1175  01=0  5+01+01 
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lire'  P7=EXP<-01.> 

1180  F'7=P7»  <:  20+M8+V8+M8+V8  >*<  20— R  >  / 100000000 

1181  GOSUB  17 
1135  GOTO  1135 
1130  P7=0 

1135  P  =P  7 + P  8 * P 3 / J 0 

1200  PRINT  10000O* I +J,  100080*K-M_j  P 

1201  PRINT  5* < M— M0 > /It0*  5*< V-V0>/V0,  J0 

1202  PRINT 

1205  P6=P6+P 

1206  IF  L>0THEN  1235 
1210  IF  M<M2THEN  1225 

1215  F'6=P6/<.  3973*.  3372  > 

1216  P1=P1+P6 
1220  F'2=P2+F'6+P6 
1220  GOTO  130  0 

1225  IF  M:>N2THEN  1265 
1236  IF  L=C3THEN  1020 
1240  IF  V< VI THEN  1020 
1245  GOTO  104O 
1265  F3=L-1  GOTO  1020 
1200  RETURN 
400O  OEFFN ' 77 

4005  PEN  CALCULATE  IF  V>V2 

4006  REM  CALCULATE  V 
4010  S1=0 

4015  51=51+1 
4020  S2=S1+V2 
4025  G  0  S  U  B  "  3  3  <  5  2  > 

4030  IF  0400 THEN  4015 
4035  GOSUB  •"34 
4040  RETURN 
5100  OEFFN  ■  10 •'  01  > 

5105  PEN  SAMPLE  FROM  NORMAL  DISTRIBUTION 

5110  GOSUB  '33 

5115  02=01+ <0  4+1.  2*Z3> 

5120  GOSUB  "93 

5125  02  =65+01+0.  9972 

5126  02=5+Z3/G3 

5120  04=5  0+ 02-01 VQ1 
5125  04*EXP<-0.  5+04 "2 > 

5140  04=5+04,-'  <0.  9972+65*01 > 

5145  IF  04 <02 THEN  5110 
5150  RETURN 
6000  DEFFN ’ 67 

6005  REM  PICK  AND  DEFINE  SHOT  LINES 

6007  J=.J  +  1 

6010  GOSUB  ' 99 

6012  S2<  1  >  =C+  2+29-1  > 

6  0 15  G  0  S  U  B  •"  9  9 


n  r  j  f  j 


6255 

GOTO  6340 

6260 

Ol=SGRv01  :■ 

*5265 

X  5=0  5  +  <  -  B9+ 0 1 > / A 9 

6270 

.76=0  5+C-B9-I 

21  :>  /  A  9 

6.27*5 

V5=S1<2)+M<2 

> * < X5-S 1 < 1 > > /  W  < 1 > 

6280 

V6=S1  <  2  >  +w  <:  2 

>  *  <  x6-si  <  i  >  >  /w  <  i  > 

6285 

Z5=S1  <  3 +H  v  3 

>  *  <  x5-si  <  i  >  :>  /w  <  l  > 

6280 

Z6=S1  <  3  >  +W  v  3 

y +  <  x6-s  i  <:  i  >  >  /w  <:  l  > 

6295 

R3=<SK1>-X5 

--2+ <  si  <  2  > -vs:.'  -2+ < 

6-700 

R4= < Sl< 1 > —76 

y  ~2+  < SIC 2 >  — V6 >  "2+ < 

62.05 

IF  R2CR4THEN 

6315 

6310 

X5=X6 :  V5=V6 

:  Z5=Z6 

6315 

R=Y5+V5+Z5+Z' 

5 

6320 

R=S0R  <  R01 

6340 

RETURN 

7 000 

DEFFN  ■'  33  <  01  > 

700S  REM  CALCULATE  F  V  > 


76*10 

02=Cl*E5+:  1  01“ 

C5> 

7015 

02=  <  01— VS  y  .'Q2 

7020 

02=02 'E4 

7025 

04=02  ~E3 

70i76* 

0  4=C6+E  6  +  C01' 

C£i  •’  +04 

76*2  5 

1714=61 3— 04 

76*401 

Q4=M8— 04 

701501 

RETURN 

7200* 

DEFFN  24 

7205 

REM  BINARY  SEARCH 

7210 

06=1.  000001+S 

2 

7212 

IF  S1MTHEN  7 

215 

7214 

05=0i  9 9999 8 

:  GOTO  7220 

7215 

05=0.  99 9 9 9 *  <  S 

1-1 > *V2 

7220 

07=0.  5* <05+06 

> 

r' 

IF  < 0 6 — 0 5 y / 0 6 

06i  0000O1THEN 

7225 

GOSUB  '22  C  07  > 

7230 

IF  04  :01.  ©THEN 

7240 

7225 

06=07 

7228 

GOTO  7226* 

72401 

Q5=07 

7245 

GOTO  72201 

7276* 

RETURN 

86*06i 

DEFFN  "  17 

86*6i5 

REM  CALCULATE 

JACOBIAN 

86*101 

01 =M “ E2 

8015 

02=E6 

80*201 

02=VC0 

8025 

04= C6 

801201 

05=0 1/M 

86*25 

06=E5 

86*46* 

07 =M”E1 

86*45 

08=01 
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U'l  l/l 


SO50 

09=Q7/M 

5055 

i7!0=V*r:5 

5060 

J  1=1— E2*Q4*Q2*Q5*0  S 

8065 

J2=-C0*Q4+Q 

2*01*02.  ‘V 

•?O70 

T2=-E1*Q8*0 

6*09*00 

3075 

J4=1-C5*Q8* 

O6*Q7*O0/V 

8080 

J  0= Jl* J4— J 2 

:*.J3 :  T0=ABSO_T0> 

8150 

RETURN 

8000 

DEFFN  89 

8005 

REM  RANDOM 

NUMBER  MAIN  SR 

8010 

25=25*28 

8015 

GOSUB  '98 

8020 

28=5*23 

8025 

GOSUB  "98 

8020 

29=28/67108 

864 

8025 

RETURN 

8050 

DEFFN  98 

8055 

REM  RANDOM 

NUMBER  SECONDER' 

30*£*0 

05=2 8 .■•'671 0  8 

864 

30*5 

06=67 108364 

*  I  NT  05  "> 

•3070 

28=28-06 

3075 

RETURN 
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